#Jupyter magic command to improve teh display of charts in the Notebook
%matplotlib inline

3. Scenario analysis

An essential feature of a model is that when given a specific set of inputs (the exogenous variables to the model) it will always return the same results.

Below a new ModelFlow sessionis prepared, initializing a pandas session and importing and solving a saved WBG model (NB: these are precisely the same commands they used to start the previous chapter) and would form the essential initialization commands of any python session using ModelFlow.

# import the model class from modelflow package
from modelclass import model 
import modelmf       # Add useful features to pandas dataframes 
                     # using utlities initially developed for modelflow

model.widescreen()   # These modelflow commands ensure that outputs from modelflow play well with Jupyter Notebook
model.scroll_off()

%load_ext autoreload   
%autoreload 2


#Load a saved version of the Pakistan model and solve it, 
#saving the results in the model object mpak, and the resulting dataframe in bline

#Replace the path below with the location of the pak.pcim file on your computer
mpak,bline = model.modelload('..\models\pak.pcim', \
                                alfa=0.7,run=1,keep= 'Baseline')

file read:  C:\modelflow manual\papers\mfbook\content\models\pak.pcim

As noted, when the model is solved without changing any inputs (as was the case of the load) the model should return (reproduce) exactly the same data as before[^fn2]. To test this for mpak the results from the simulation can be compared by using the basedf and lastdf DataFrames.

[^fn2:] If it does not, the model has violated he principle of reproducibility and there is something wrong (usually one of the identities does not hold).

Below, the percent difference between the values of the variables for real GDP and Consumer demand in the two dataframes .basedf and lastdf is zero following a simulation where the inputs were not changed – confirming the reproduction of results.

# Need statement to change the default format
mpak.smpl(2020,2030)
mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN'].difpctlevel.mul100.df

	PAKNYGDPMKTPKN	PAKNECONPRVTKN
2020	0.0	0.0
2021	0.0	0.0
2022	0.0	0.0
2023	0.0	0.0
2024	0.0	0.0
2025	0.0	0.0
2026	0.0	0.0
2027	0.0	0.0
2028	0.0	0.0
2029	0.0	0.0
2030	0.0	0.0

3.1. Different kinds of simulations

The modelflow package performs 4 different kinds of simulation:

1. A shock to an exogenous variable in the model

2. An exogenous shock of a behavioural variable, executed by exogenizing the variable

3. An endogenous shock of a behavioural variable, executed by shocking the add factor of the variable.

4. A mixed shock of a behavioural variable, achieved by temporarily exogenixing the variable.

Although technically modelflow would allow us to shock identities, that would violate their nature as accounting rules. Effectively such a shock would break the economic sense of the model.

As a result, this possibility is not discussed.

3.1.1. A shock to an exogenous variable

A World Bank model will reproduce the same values if inputs (exogenous variables) are not changed. In the simulation below, the oil price is changed – increasing by $25 for the three years between 2025 and 2027 inclusive.

3.1.1.1. Preparing the data for simulation

The following steps are performed to prepare the simulation:

1. A new input dataframe is created as a copy of the original.

2. The oil price in the new dataframe is modified using the mfcalc method for the three years in question.

3. Finally pandas math is used to display the initial value, the changed value and the difference between the two, confirming that the mfcalc statement revised the oil price data.

#Make a copy of the baseline dataframe
oilshockdf=mpak.basedf
oilshockdf=oilshockdf.mfcalc("<2025 2027> WLDFCRUDE_PETRO = WLDFCRUDE_PETRO +25")

compdf=mpak.basedf.loc[2000:2030,['WLDFCRUDE_PETRO']]
compdf['LASTDF']=oilshockdf.loc[2000:2030,['WLDFCRUDE_PETRO']]
compdf['Dif']=compdf['LASTDF']-compdf['WLDFCRUDE_PETRO']

compdf.loc[2024:2030]

	WLDFCRUDE_PETRO	LASTDF	Dif
2024	80.367180	80.367180	0.0
2025	85.336809	110.336809	25.0
2026	90.613742	115.613742	25.0
2027	96.216983	121.216983	25.0
2028	102.166709	102.166709	0.0
2029	108.484346	108.484346	0.0
2030	115.192643	115.192643	0.0

3.1.1.2. Running the simulation

Having created a new dataframe comprised of all the old data plus the changed data for the oil price, a simulation can now be run.

In the command below, the simulation is run from 2020 to 2040, using the oilshockdf as the input DataFrame. The results of the simulation are assigned to a new DataFrame named ExogOilSimul. The Keep command ensures that the mpak model object stores (keeps) a copy of the results identified by the text name ‘$25 increase in oil prices 2025-27’.

#Simulate the model 
ExogOilSimul = mpak(oilshockdf,2020,2040,keep='$25 increase in oil prices 2025-27') 

3.1.1.2.1. Results

ModelFlow tools can be used to visualize the impacts of the shock; as a print out; as charts and within Jupyter notebook as an interactive widget.

The display below confirms that the shock was executed as desired. The dif.df method returns the difference between the .lastdf and .basedf values of the selected variable(s) as a DataFrame. The with mpak.set_smpl(2020,2030): clause temporarily restricts the sample period over which the following indented commands are executed.

Alternatively the mpak.smpl(2020,2030)could be used. This would restricts the time period of over which all subsequent commands are executed.

with mpak.set_smpl(2020,2030):
    print(mpak['WLDFCRUDE_PETRO'].dif.df);

      WLDFCRUDE_PETRO
2020              0.0
2021              0.0
2022              0.0
2023              0.0
2024              0.0
2025             25.0
2026             25.0
2027             25.0
2028              0.0
2029              0.0
2030              0.0

Below the impact of this change on a few variables are expressed graphically and in a table.

The first variable PAKNYGDPMKTPKN is Pakistan’s real GDP, the second PAKNECONPRVTKN is real consumption and the third is the Consumer price deflator PAKNECONPRVTXN.

The modifier .difpctlevel.mul100 instructs modelflow to calculate the difference in the selected variables and express them as a percent of the original level.

\[ = \frac{x_{new}-x_{original}}{x_{original}} * 100 \]

mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="Impact of temporary $25 hike in oil prices")

with mpak.set_smpl(2020,2035):
    print(round(mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.df,2));

      PAKNYGDPMKTPKN  PAKNECONPRVTKN  PAKNEIMPGNFSKN  PAKNECONPRVTXN
2020            0.00            0.00            0.00            0.00
2021            0.00            0.00            0.00            0.00
2022            0.00            0.00            0.00            0.00
2023            0.00            0.00            0.00            0.00
2024            0.00            0.00            0.00            0.00
2025           -0.89           -1.32           -1.49            1.64
2026           -0.85           -1.48           -2.65            1.35
2027           -0.64           -1.37           -3.19            1.08
2028            0.34           -0.08           -2.17           -0.51
2029            0.50            0.20           -1.25           -0.43
2030            0.45            0.19           -0.80           -0.31
2031           -0.04            0.02           -0.10            0.26
2032           -0.06            0.01           -0.01            0.20
2033           -0.08            0.00            0.03            0.15
2034           -0.08            0.01            0.04            0.11
2035           -0.07            0.03            0.05            0.08

The graphs show the change in the level as a percent of the previous level. They suggest that a temporary $25 oil price hike would reduce GDP in the first year by about 0.9 percent, that the impact would diminish by the third year to -.64 percent, and then turn positive in the fourth year when the price effect was eliminated.

The impacts on household consumption are stronger but follow a similar pattern.

The GDP impact is smaller partly because the decline in domestic demand reduces imports. Because imports enter into the GDP identity with a negative sign, a reduction in imports actually increase aggregate GDP – or in this case partially offsets the declines coming from reduced consumption (and investment).

Finally as could be expected, initially prices rise sharply with higher oil prices. However, as the slow down in growth is felt, inflationary pressures turn negative and the overall impact on the price level turns negative. The graph and table above shows what is happening to the price level. To see the impact on inflation (the rate of growth of prices), a separate graph can be generated using difpct.mul100, which shows the change in the rate of growth of variables where the growth rate is expressed as a per cent \(\bigg[\bigg(\frac{x^{shock}_t}{x^{shock}_{t-1}}-1\bigg)\) \( - \bigg(\frac{x^{baseline}_t}{x^{baseline}_{t-1}}-1\bigg)\Bigg]*100\).

mpak['PAKNECONPRVTXN'].difpct.mul100.plot(title="Change in inflation from a temporary $25 hike in oil prices")

Ib how come this graph shows up so small. How can we affect its size?

This view, gives a more nuanced result. The inflation rate increases initially by about 1.2 percentage points, but falls compared with the baseline below in the 2026-2027 period as the as the influence of the slowdown in GDP more than offsets the continued inflationary impetus from the lagged increase in oil prices. In 2028, when oil prices drop back to their previous level, there is an additional dis-inflationary force and sharp drop in inflation as compared with the baseline. Overtime, the boost to demand from lower prices translates into an acceleration in growth and a return of inflation back to its trend rate.

3.1.2. An exogenous shock to a Behavioural variable

Behavioural equations can be de-activated by exogenizing them, either for the entire simulation period, or for a selected sub-period. In this example, consumption is exogenized for the entire simulation period.

To motivate the simulation, it is assumed that a change in weather patterns has increased the number of sunny days by 10 percent. This increases households happiness and causes them to permanently increase their spending by 2.5% beginning in 2025.

Such a shock can be specified either manually or by using the.fix() method. Below the simpler .fix() method is used, but the equivalent manual steps performed by .fix() are also explained.

To exogenize PAKNECONPRVTKN for the entire simulation period, initially a new DataFrame Cfixed is created as a slightly modified version of mpak.basedf using the .fix() command.

Cfixed=mpak.fix(mpak.basedf,PAKNECONPRVTKN)

This does two things, that could have been done manually. First it sets the dummy variable PAKNECONPRVTKN_D=1 for the entire simulation period. Recall the consumption equation like all behavioural equations of World Bank models implemented in ModelFlowis expressed in two parts.

\[ cons= (1-cons_D)*\bigg[C'(X)\bigg] + cons_d*cons_x\]

When \(cons_D=1\) (as it does in this scenario) the first part of the equation \((1-cons_D)*C'(X)\) evaluates to zero and consumption is equal to (1)* \(cons_x\). If instead (which would be the normal case \(cons_d\) were set to zero. the the equation would simplify to \( cons= C'(X) \)

Then .fix() method then sets the variable PAKNECONPRVTKN_X in the Cfixed dataframe equal to the value of PAKNECONPRVTKN in the basedf .DataFrame. All the other variables are just copies of their values in .basedf.

With PAKNECONPRVTKN_D=1 throughout, the normal behavioral equation is effectively de-activated or exogenized and consumption will just equal its exogenized value: \(PAKNECONPRVTKN=PAKNECONPRVTKN_X\).

mpak.smpl() # reset the active sample period to the full model.
Cfixed=mpak.fix(bline,'PAKNECONPRVTKN')

The folowing variables are fixed
PAKNECONPRVTKN

For the moment, the equation is exogenized but the values have been set to the same values as the .basedf dataframe, so solving the model will not change anything.

The .upd() method can now be used to implement the assumption that Real consumption ( PAKNECONPRVTYKN) would be 2.5% stronger.

Cfixed=Cfixed.upd("<2025 2040> PAKNECONPRVTKN_X  * 1.025")

To perform the simulation, the revised CFixed DataFrame is input to the mpak model solve routine.

CFixedRes = mpak(Cfixed,2020,2040,keep='2.5% increase in C 2025-40 (fix)')

And then the results can be examined graphically as before.

CFixedRes = mpak(Cfixed,2020,2040,keep='2.5% increase in C 2025-40') # simulates the model 
mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="Impact of a permanent 2.5% increase in Consumption")

Below results are displayed in tabular form. Note the use of the pandas options with the with clause. This sets the display format of floating point variables to two decimal points. The second with clause temporarily restricts the display to the period 2020 to 2040.

import pandas as pd
with pd.option_context('display.float_format', '{:,.2f}'.format):
    with mpak.set_smpl(2020,2040):
        print(mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.df)

      PAKNYGDPMKTPKN  PAKNECONPRVTKN  PAKNEIMPGNFSKN  PAKNECONPRVTXN
2020            0.00            0.00            0.00            0.00
2021            0.00            0.00            0.00            0.00
2022            0.00            0.00            0.00            0.00
2023            0.00            0.00            0.00            0.00
2024            0.00            0.00            0.00            0.00
2025            2.01            2.50            2.27            0.44
2026            2.07            2.50            2.43            1.06
2027            2.05            2.50            2.59            1.69
2028            1.99            2.50            2.78            2.31
2029            1.92            2.50            2.99            2.90
2030            1.83            2.50            3.22            3.47
2031            1.43            2.50            4.03            4.53
2032            1.37            2.50            4.18            4.92
2033            1.30            2.50            4.34            5.29
2034            1.23            2.50            4.50            5.64
2035            1.16            2.50            4.66            5.97
2036            1.09            2.50            4.81            6.28
2037            1.03            2.50            4.96            6.56
2038            0.96            2.50            5.10            6.82
2039            0.90            2.50            5.24            7.06
2040            0.84            2.50            5.36            7.28

The permanent rise in consumption by 2.5 percent causes a temporary increase in GDP of close to 2% (1.86). Higher imports tend to diminish the effect on GDP. Over time higher prices due to the inflationary pressures caused by the additional demand cause the GDP impact to diminish to close to less than 1 percent by 2040.

3.1.3. Exogenize a behavioural variable and temporarily shock it

The third method of formulating a scenario involves exogenizing an endogenous variable and shocking its value for a defined period of time. The methodology is the same except the period for which the variable is exogenized is different.

Here the set up is basically the same as before.

#reset the active sample period to the full period
mpak.smpl()                                  
# create a copy of the bline DataFrame, but setting the PAKNECONPRVTKN_D variable to 1 for the period 2025 through 2027
CTempExogAll=mpak.fix(bline,'PAKNECONPRVTKN') 
# multiply the exogenized value of consumption by 2.5% for 2025 through 2027
CTempExogAll=CTempExogAll.upd("<2025 2027> PAKNECONPRVTKN_X * 1.025")  

#Solve the model
CTempXAllRes = mpak(CTempExogAll,2020,2040,keep='2.5% increase in C 2025-27 -- exog whole period') # simulates the model 
mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="Temporary hike in Consumption 2025-2027")

The folowing variables are fixed
PAKNECONPRVTKN

The results are quite different. GDP is boosted initially as before but when consumption drops back to its pre-shock level, GDP and imports decline sharply.

Prices (and inflation) are higher initially but when the economy starts to slow after 2025 prices actually fall (deflation). While prices are falling, the level of prices remains higher at the end of the simulation.

3.1.3.1. Temporary shock exogenized for the whole period

This scenario is the same as the previous, but this time the --KG (keep_growth) option is used to maintain the pre-shock growth rates of consumption in the post-shock period. Effectively this is the same as a permanent increase in the level of consumption by 2.5% because the final shocked value of consumption (which was 2.5% higher then its pre-shock level) is grown at the same pre-shock rate – ensuring that all post-shock variables are also up by 2.5%.

mpak.smpl() # reset the active sample period to the full model.
CTempExogAllKG=mpak.fix(bline,'PAKNECONPRVTKN')
CTempExogAllKG = CTempExogAllKG.upd('''
<2025 2027> PAKNECONPRVTKN_X * 1.025 --kg
''',lprint=0)

#Now we solve the model
CTempXAllResKG = mpak(CTempExogAllKG,2020,2040,keep='2.5% increase in C 2025-27 -- exog whole period --KG=True') # simulates the model 
mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="2.5% boost to cons 2025-27 --kg=True")

The folowing variables are fixed
PAKNECONPRVTKN

3.1.4. Exogenize the variable only for the period during which it is shocked

This scenario introduces a subtle but import difference. Here the variable is again exogenized using the fix syntax. But this time it is exogonized only for the period where the variable is shocked.

This means that the consumption equation will only be de-activated for three years (instead of the whole period as in the previous examples). As a result, the values that consumption takes in 2028, 2029, … 2040 depend on the model, not the level it was set to when exogenized (which was the case in the 3 previous versions).

Looking at the maths of the model the consumption equation is effectively split into two.

for the period before 2025 \(cons_d=0\) and the consumption equation simplifies to:

\(cons=C(X)\)

for the period 2025-2028 it is exogenized (\(cons_d=1\)) so it simplifies to:

\(cons=cons_x\)

but in the final period 2028-2040 (\(cons_d=0\)) and the equation reverts to:

\(cons=C(X)\)

mpak.smpl() # reset the active sample period to the full model.
CExogTemp=mpak.fix(bline,'PAKNECONPRVTKN',2025,2027)                             #Consumption is exogenized only for three years 2025 2026 and 2027 PAKNECONPRVTKN_D=1 for 2025,2026, 2027 0 elsewhere.
CExogTemp = CExogTemp.upd('<2025 2027> PAKNECONPRVTKN_X * 1.025',lprint=0)       #In subsequent years it's level will be determined by the equation 

#Solve the model
CExogTempRes = mpak(CExogTemp,2020,2040,keep='2.5% increase in C 2025-27 -- temporarily exogenized') # simulates the model 
mpak['PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="Temporary 2.5% boost to cons 2025-27 - equation active")

The folowing variables are fixed
PAKNECONPRVTKN

These results have subtle differences compared with the previous. The most obvious is visible in looking at the graph for Consumption. Rather than reverting immediately to its earlier pre-shock level, it falls more gradually and actually overshoots (falls below its earlier level), before returning slowly to its pre-shock level. That is because unlike in the previous shocks, its path is being determined endogenously and reacting to changes elsewhere in the model, notably changes to prices, wages and government spending as well as the lagged level of consumption.

print('Consumption base and shock levels\r\n');

print('Real values in 2030');
print(f'Base value:  {bline.loc[2028,"PAKNECONPRVTKN"]:,.0f}.\tShocked value: {CExogTempRes.loc[2028,"PAKNECONPRVTKN"]:,.0f}.\r\n'
    f'Percent difference: {round(100*((CExogTempRes.loc[2030,"PAKNECONPRVTKN"]-bline.loc[2028,"PAKNECONPRVTKN"])/bline.loc[2028,"PAKNECONPRVTKN"]),2)}')
print('\r\nReal values in 2040');
print(f'Base value:  {bline.loc[2040,"PAKNECONPRVTKN"]:,.0f}.\tShocked value: {CExogTempRes.loc[2040,"PAKNECONPRVTKN"]:,.0f}.\r\n'
    f'Percent difference: {round(100*((CExogTempRes.loc[2040,"PAKNECONPRVTKN"]-bline.loc[2040,"PAKNECONPRVTKN"])/bline.loc[2040,"PAKNECONPRVTKN"]),2)}')

Consumption base and shock levels

Real values in 2030
Base value:  27,241,278.	Shocked value: 27,616,949.
Percent difference: 5.36

Real values in 2040
Base value:  38,676,995.	Shocked value: 38,693,167.
Percent difference: 0.04

3.1.5. Simulation with Add factors

Add factors are a crucial element of the macromodels of the World Bank and serve multiple purposes.

In simulation, add-factors allow simulations to be conducted without de-activating behavioural equations. Such shocks are often referred to as endogenous shocks because the equation of the behavioural variable that is shocked remains active throughout.

In some ways they are very similar to a temporary exogenous shock. Both ways of producing the shock allow the shocked variable to respond endogenously in the period after the shock. The main difference between the two approaches is that:

• Endogenous shocks (Add-Factor shocks) allow the shocked variable to respond to changed circumstances that occur during the period of the shock.

  • This approach makes most sense for “animal spirits”, shocks where the underlying behaviour is expected to change.

  • It also makes sense when actions of one part of an aggregate is likely to impact behaviour of other sectors within an aggregate

  • increased investment by a particular sector would be an example here as the associated increase in activity is likely to increase investment incentives in other sectors, while increased demand for savings will increase interest rates and the cost of capital operating in the opposite direction.

  • Sustained changes in behaviour, for example increased propensity to invest because of improved recognition

• Exogenous shocks to endogenous variables fix the level of the shocked variable during the shock period.

  • Changes in government spending policy, something that is often largely an economically exogenous decision.

3.1.5.1. Simulating the impact of a planned investment

The below simulation uses the add-factor to simulate the impact of a 3 year investment program beginning in 2025 of 1 percent of GDP per year, that is financed through an increase in foreign direct investment. This might reflect a specific large scale plant that is being constructed due to a deal reached by the government with a foreign manufacturer. The add-factor approach is chosen because the additional investment is likely to increase demand for the products of other firms, which is likely to incite them to add to their investments as well – both after the shock as in previous examples but also during the period that investment is being shocked.

3.1.5.1.1. How to translate the economic shock into a model shock

Add-factors in the MFMod framework are applied to the intercept of an equation (not the level of the dependent variable). This preserves the estimated elasticities of the equation, but makes introduction of an add-factor shock somewhat more complicated than the exogenous approach. Below a step-by-step how-to guide:

1. Identify numerical size of the shock

2. Examine the functional form of the equation, to determine the nature of the add factor. If the equation is expressed as a:

   • growth rate then the add-factor will be an addition or subtraction to the growth rate

   • percent of GDP (or some other level) then the add-factor will be an addition or subtraction to the share of growth.

   • Level then the add-factor will be a direct addition to the level of the dependent variable

3. Convert the economic shock into the units of the add-factor

4. Shock the add-factor by the above amount and run the model

   • Note the add-factor is an exogenous variable in the model, so shocking it follows the well established process for shocking an exogenous variable.

3.1.5.1.2. Determine the size of shock

Above the shock was identified as a 1 percent of GDP increase in FDI that flows directly into private-sector investment. A first step would be to determine the variables that need to be shocked (FDI) and private investment. To do this we can query the variable dictionary.

mpak['!*invest*'].des

PAKNEGDIFGOVKN        : Pub investment real
PAKNEGDIFPRVKN        : Prvt. Investment real
PAKNEGDIFPRVKN_A      : Add factor:Prvt. Investment real
PAKNEGDIFPRVKN_D      : Fix dummy:Prvt. Investment real
PAKNEGDIFPRVKN_FITTED : Fitted  value:Prvt. Investment real
PAKNEGDIFPRVKN_X      : Fix value:Prvt. Investment real
PAKNEGDIFTOTKN        : Investment real

Querying for all variables that “invest” in their descriptor gives us the private investment variable.

3.1.5.1.3. Identify the functional form(s)

To understand how to shock using the add factor, it is essential to understand how the add-factor enters into the equation.

Addfactor is on the intercept of	Shock needs to be calculated as
a growth equation	a change in the growth rate
Share of GDP	a percent of GDP
Level	as change in the level

Use the .eviews command command to identify the functional forms of the equation to be shocked.

mpak['PAKNEGDIFPRVKN'].eviews

PAKNEGDIFPRVKN : (PAKNEGDIFPRVKN/PAKNEGDIKSTKKN( - 1)) = 0.00212272413966296 + 0.970234989019907*(PAKNEGDIFPRVKN( - 1)/PAKNEGDIKSTKKN( - 2)) + (1 - 0.970234989019907)*(DLOG(PAKNYGDPPOTLKN) + PAKDEPR) + 0.0525240494260597*DLOG(PAKNEKRTTOTLCN/PAKNYGDPFCSTXN)

In this case the equation is written as a share of the capital stock in the preceding period PAKNEGDIKSTKKN(-1).

3.1.5.1.4. Calculate the size of the required add factor shock

The shock to be executed is 0.5 percent of GDP.

It is assumed that the financing will come from FDI and that all the money will be spent in one year on private investment.

The private investment equation is written as a share of the capital stock. Therefore, the add-factor needs to be shocked by adding 1 percent of GDP to private investment in 2028 divided by the capital stock in 2028.

#Create a DataFrame AFShock that is equal tothe baseline
AFShock=bline

#Display the level of the AF
print("Pre shock levels")
AFShock.loc[2025:2030,['PAKNEGDIFPRVKN_A','PAKNEGDIFPRVKN','PAKNEGDIKSTKKN']]

#print(AFShock.loc[2025:2030,'PAKNEGDIFPRVKN']/AFShock.loc[2025:2030,'PAKNYGDPMKTPKN']*100)

Pre shock levels

	PAKNEGDIFPRVKN_A	PAKNEGDIFPRVKN	PAKNEGDIKSTKKN
2025	-0.000458	1.602854e+06	4.730392e+07
2026	-0.000389	1.581104e+06	4.814879e+07
2027	-0.000331	1.569541e+06	4.900980e+07
2028	-0.000281	1.569141e+06	4.989869e+07
2029	-0.000239	1.580577e+06	5.082694e+07
2030	-0.000203	1.604394e+06	5.180590e+07

Below the mfcalc routine is used to set the addfactor variable equal to its previous value plus the equivalent of 1 percent of GDP when expressed as a percent of the previous period’s level of capital stock.

AFShock=AFShock.mfcalc("<2028 2028> PAKNEGDIFPRVKN_A = PAKNEGDIFPRVKN_A + (.01*PAKNYGDPMKTPKN/PAKNEGDIKSTKKN(-1))");

print("Shocked AF levels")
AFShock.loc[2025:2030,'PAKNEGDIFPRVKN_A']

Shocked AF levels

2025   -0.000458
2026   -0.000389
2027   -0.000331
2028    0.005884
2029   -0.000239
2030   -0.000203
Name: PAKNEGDIFPRVKN_A, dtype: float64

3.1.5.1.5. Run the shock

AFShockRes = mpak(AFShock,2020,2040,keep='1% of GDP increase in FDI and private investment (AF shock)')
mpak['PAKNYGDPMKTPKN PAKNEGDIFPRVKN PAKNECONPRVTKN PAKNEIMPGNFSKN PAKNEGDIFTOTKN PAKNECONPRVTXN'].difpctlevel.mul100.plot(title="Add factor shock on private investment 1% of GDP")

The above graphs are expressed as a percent of the baseline value of each value. As a result, the shock of 1% of GDP, expressed as a percent of private investment (which is much smaller than GDP) appears quite large.

Below, to double-check two variables IFPRVOLD_ORIG and IFTOT_ORIG are created that reflect pre-shock private and total investment as a share of pre-shock GDP. Two additional variables IFPRV_SHOCK and IFTOT_SHOCK are also created as the shocked values of private and total investment as a percent of the original GDP. The difference between the tqo sets of variables is the increase in fixed private and fixed total investment as a per cent of the same denominator (the original level of GDP).

The ex post change in private investment and of total investment in 2028 is 1.06 and 1.11 percent respectively, 0.6 and 1.1 percentage points larger than the actual shock 1 percent of GDP shock to the addfactor.

This difference represents the endogenous reaction of other investors in the same time period to the changed circumstances. It is precisely to capture this effect that an endogenous or add-factor shock is employed.

AFShockRes['GDPOLD']=bline['PAKNYGDPMKTPKN']
AFShockRes['IFPRVOLD']=bline['PAKNEGDIFPRVKN']
AFShockRes['IFTOTOLD']=bline['PAKNEGDIFTOTKN']
AFShockRes=AFShockRes.mfcalc('''
                       IFPRV_SHOCK = PAKNEGDIFPRVKN / GDPOLD*100
                       IFTOT_SHOCK = PAKNEGDIFTOTKN / GDPOLD*100
                       IFPRV_Orig = IFPRVOLD / GDPOLD*100
                       IFTOT_Orig = IFTOTOLD / GDPOLD*100
                       ''')

print(round(AFShockRes.loc[2020:2030,['PAKNEGDIFPRVKN_','IFPRV_OLD','PAKNEGDIFTOTKN_','IFTOT_OLD']],2))

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[23], line 11
      3 AFShockRes['IFTOTOLD']=bline['PAKNEGDIFTOTKN']
      4 AFShockRes=AFShockRes.mfcalc('''
      5                        IFPRV_SHOCK = PAKNEGDIFPRVKN / GDPOLD*100
      6                        IFTOT_SHOCK = PAKNEGDIFTOTKN / GDPOLD*100
      7                        IFPRV_Orig = IFPRVOLD / GDPOLD*100
      8                        IFTOT_Orig = IFTOTOLD / GDPOLD*100
      9                        ''')
---> 11 print(round(AFShockRes.loc[2020:2030,['PAKNEGDIFPRVKN_','IFPRV_OLD','PAKNEGDIFTOTKN_','IFTOT_OLD']],2))

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:1097, in _LocationIndexer.__getitem__(self, key)
   1095     if self._is_scalar_access(key):
   1096         return self.obj._get_value(*key, takeable=self._takeable)
-> 1097     return self._getitem_tuple(key)
   1098 else:
   1099     # we by definition only have the 0th axis
   1100     axis = self.axis or 0

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:1289, in _LocIndexer._getitem_tuple(self, tup)
   1286 if self._multi_take_opportunity(tup):
   1287     return self._multi_take(tup)
-> 1289 return self._getitem_tuple_same_dim(tup)

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:955, in _LocationIndexer._getitem_tuple_same_dim(self, tup)
    952 if com.is_null_slice(key):
    953     continue
--> 955 retval = getattr(retval, self.name)._getitem_axis(key, axis=i)
    956 # We should never have retval.ndim < self.ndim, as that should
    957 #  be handled by the _getitem_lowerdim call above.
    958 assert retval.ndim == self.ndim

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:1332, in _LocIndexer._getitem_axis(self, key, axis)
   1329     if hasattr(key, "ndim") and key.ndim > 1:
   1330         raise ValueError("Cannot index with multidimensional key")
-> 1332     return self._getitem_iterable(key, axis=axis)
   1334 # nested tuple slicing
   1335 if is_nested_tuple(key, labels):

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:1272, in _LocIndexer._getitem_iterable(self, key, axis)
   1269 self._validate_key(key, axis)
   1271 # A collection of keys
-> 1272 keyarr, indexer = self._get_listlike_indexer(key, axis)
   1273 return self.obj._reindex_with_indexers(
   1274     {axis: [keyarr, indexer]}, copy=True, allow_dups=True
   1275 )

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexing.py:1462, in _LocIndexer._get_listlike_indexer(self, key, axis)
   1459 ax = self.obj._get_axis(axis)
   1460 axis_name = self.obj._get_axis_name(axis)
-> 1462 keyarr, indexer = ax._get_indexer_strict(key, axis_name)
   1464 return keyarr, indexer

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexes\base.py:5876, in Index._get_indexer_strict(self, key, axis_name)
   5873 else:
   5874     keyarr, indexer, new_indexer = self._reindex_non_unique(keyarr)
-> 5876 self._raise_if_missing(keyarr, indexer, axis_name)
   5878 keyarr = self.take(indexer)
   5879 if isinstance(key, Index):
   5880     # GH 42790 - Preserve name from an Index

File ~\miniconda3\envs\mfbook_310\lib\site-packages\pandas\core\indexes\base.py:5935, in Index._raise_if_missing(self, key, indexer, axis_name)
   5933     if use_interval_msg:
   5934         key = list(key)
-> 5935     raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   5937 not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())
   5938 raise KeyError(f"{not_found} not in index")

KeyError: "None of [Index(['PAKNEGDIFPRVKN_', 'IFPRV_OLD', 'PAKNEGDIFTOTKN_', 'IFTOT_OLD'], dtype='object')] are in the [columns]"

4. Report writing and scenario results

ModelFlow, standard pandas routines and other python libraries like Matplotlib and Plotly can be used to visualize and compare dataframes and therefore the results from scenarios – as indeed has been done in the preceding paragraphs.

In addition, ModelFlow also provides several specific routines that make such comparisons easier.

4.1. The Keep option

The Keep option facilitates the comparison of results from different scenarios run on a give model object. In each of the simulations executed above, the keep option was activated. This causes the results from each simulation to be stored in the model object in a DataFrame that can be identified by the descriptor given to it.

4.2. The .keep_plot() method

The keep_plot method can be used to plot and compare results from the various scenarios that had been run earlier using the keep= option.

By default the results across all scenarios for each selected variables will be shown on one chart at a time.

mpak.keep_plot('PAKNYGDPMKTPCN PAKNECONPRVTKN PAKNEIMPGNFSKN',  legend=True)
#show for each variable on a separate chart the results from each kept scenario

{'PAKNYGDPMKTPCN': <Figure size 1000x600 with 1 Axes>,
 'PAKNECONPRVTKN': <Figure size 1000x600 with 1 Axes>,
 'PAKNEIMPGNFSKN': <Figure size 1000x600 with 1 Axes>}

4.2.1. keep_plot() options

The variables to be displayed are listed as first argument. Variable names can include wildcards (using * for any string and ? for any character).

Transformation of data displayed:

showtype=	Use this operator
‘level’ (default)	No transformation
‘growth’	The growth rate in percent
‘change’	The yearly change (\(\Delta\))

legend placement

legend=	Use this operator
False (default)	The legends are placed at the end of the corresponding line
True	The legends are places in a legend box

Often it is useful to compare the scenario results with the baseline result. This is done with the diff argument.

diff=	Use this operator
False (default)	All entries in the keep_solution dictionary are displayed
True	The difference to the first entry is shown.

It can also be useful to compare the scenario results with the baseline result measured in percent. This is done with the diffpct argument.

diffpct=	Use this operator
False (default)	All entries in the keep_solution dictionary is displayed
True	The difference in percent to the first entry is shown

Note

'keep_plot() and .keep_plot_multi() return a python object that points to the in memory version of the rendered figure(s). This object can be used to modify the graph (see examples towards the end of this chapter.

savefig='[path/]<prefix>.<extension>' Will create a number of files with the charts.
The files will be saved location with name <path>/<prefix><variable name>.<extension>
The extension determines the format of the saved file: pdf, svg and png are the most common extensions.

4.2.2. An example using the diff=TRUE option

When diff=True (or 1) results will be shown all of the selected scenarios presented as the change in selected variables with respect to the first scenario – in this instance the scenario saved with the name baseline.

Note in this instance `baseline` and `basedf` are the same because they were defined that way.  However, there is nothing in the system that guarantees that the first `keep` scenario will be the baseline or the `basedf` scenario.

mpak.keep_plot('PAKNYGDPMKTPCN PAKNECONPRVTKN PAKNEIMPGNFSKN', diff=1, legend=True)

{'PAKNYGDPMKTPCN': <Figure size 1000x600 with 1 Axes>,
 'PAKNECONPRVTKN': <Figure size 1000x600 with 1 Axes>,
 'PAKNEIMPGNFSKN': <Figure size 1000x600 with 1 Axes>}

```{index} single: Scenarios: Keep_plot examples - showtype=growth

4.2.3. The showtype option

In this example the difference with respect first 'keep scenario baseline values are once again shown. This time the showtype option has been set to growth. As a result the data is displayed as the difference in the growth rate.

mpak.keep_plot('PAKNYGDPMKTPCN PAKNEIMPGNFSKN', diff=1,showtype='growth', legend=True)

{'PAKNYGDPMKTPCN': <Figure size 1000x600 with 1 Axes>,
 'PAKNEIMPGNFSKN': <Figure size 1000x600 with 1 Axes>}

```{index} single: Scenarios: Keep_plot examples - showtype=diffpct

4.2.4. The diffpct option

Setting diffpct=True instructs .keep_plot() to display the data as a percent deviation from the first kept scenario, which in this case was the baseline scenario.

mpak.keep_plot('PAKNYGDPMKTPCN  PAKNEIMPGNFSKN', diffpct=1,legend="Change in level as a % of first keep scenario")

{'PAKNYGDPMKTPCN': <Figure size 1000x600 with 1 Axes>,
 'PAKNEIMPGNFSKN': <Figure size 1000x600 with 1 Axes>}

4.2.4.1. Differences in percent of baseline values

In this plot, the same results are presented, but as percent deviations from the baseline values of the displayed data.

mpak.keep_plot('PAKNYGDPMKTPCN PAKNECONPRVTKN ', diffpct=1,showtype='level', legend=True)

{'PAKNYGDPMKTPCN': <Figure size 1000x600 with 1 Axes>,
 'PAKNECONPRVTKN': <Figure size 1000x600 with 1 Axes>}

4.2.5. The .keep_switch() method

The .keep_switch() method restricts the number of scenarios on which subsequent calls to .keep_plot() (and .keep_plot_multi()) are executed on. .keep_switch() can be passed a list of scenarios or using a wildcard selector.

4.2.5.1. The .keep_solutions.keys() method

The .keep_solutions.keys() method generates a list of the solutions that have been kept previously.

mpak.keep_solutions.keys()

dict_keys(['Baseline', '$25 increase in oil prices 2025-27', '2.5% increase in C 2025-40', '2.5% increase in C 2025-27 -- exog whole period', '2.5% increase in C 2025-27 -- exog whole period --KG=True', '2.5% increase in C 2025-27 -- temporarily exogenized', '1% of GDP increase in FDI and private investment (AF shock)'])

To specify exactly which scenarios to show in a keep_plot, the scenarios= option of .keepswitch() must be initialized with a “|” delimited string of the names of the scenarios (retrieved above) that are to be displayed.

By placing the .keepswitch() in a with clause the scenario restriction will only apply to indented lines that follow the with construct.

with mpak.keepswitch(scenarios='2.5% increase in C 2025-40|2.5% increase in C 2025-27 -- exog whole period|2.5% increase in C 2025-27 -- exog whole period --KG=True|2.5% increase in C 2025-27 -- temporarily exogenized'):
    mpak.keep_plot('PAKNYGDPMKTPKN PAKGGBALOVRLCN PAKGGDEBTTOTLCN',diff=False,showtype='growth',legend=True);

4.2.5.2. Keepswitch with wildcard selection

Below a series of plots are generated using a wildcard selector in the keepswitch clause.

with mpak.keepswitch(scenarios='*2025*'):
    mpak.keep_plot('PAKNYGDPMKTPKN PAKNECONPRVTKN',showtype='growth',legend=True);

4.3. The .keep_plot_multi() method

The .keep_plot_multi() method allows several charts to be displayed in a grid. The size of each chart can be set with the size=(w,h) option, where the units of with and height are in centimetres.

with mpak.set_smpl(2000,2040):
    with mpak.keepswitch(scenarios="baseline *exog*"):
        var_figs = mpak.keep_plot_multi('PAKNYGDPMKTPKN PAKNECONPRVTKN PAKNEGDIFTOTKN PAKNEIMPGNFSKN PAKNEEXPGNFSKN',2010,2100,keep_dim=0,legend=1
                                ,size=(20,20) ,diffpct=True,title='Scenario Comparison'  );

As indicated earlier both keep_plot() and keep_plot_multi() return a variable that can be used to embellish or modify the figures produced by the automatic routines.

For example the charts can be resized.

var_figs.set_size_inches(15,10)
var_figs

Individual charts can be deleted from the grid.

Note

The grid representation of the individual charts is returned as a 0-based vector of charts. Thus the first figure is the zeroeth and the second is the first.

4.3.1. Delete a chart from the grid

A chart can be deleted from the grid by referencing it and calling the .remove() method.

var_figs.axes[1].remove()
var_figs

The same mechanism can be used to revise the titles of the indidividual charts and annotate them.

var_figs.axes[0].set_title('Impact of a 2.5% increase in C using exog method');    # many properties can be set afterward 
var_figs.axes[1].set_title('Impact of a 2.5% increase in C using temporary exog method');
var_figs

var_figs.axes[0].set_xlabel('Year')
var_figs.axes[0].set_ylabel('Change percent\nchange in level',fontsize=10)
var_figs.axes[0].yaxis.set_label_coords(-0.1,1.02)
var_figs.axes[0].xaxis.set_label_coords(.95,-.06)
var_figs

Variables pointing to the individual charts can be defined and used to make modifications to individual charts within the overall figure.

fig1=var_figs.axes[0]
fig2=var_figs.axes[1]

fig2.set_ylabel('Percent change\nin level',fontsize=10)
fig2.yaxis.set_label_coords(-0.1,1.02) #place axes labels
fig2.xaxis.set_label_coords(.95,-.06)



fig2.text(2040.,0.4, 'Some text in a box', 
          color='yellow',bbox=dict(facecolor='red', alpha=0.5));
fig2.text(2040.,0.1, 'Some nice text', 
          style='italic',color='green');
          
var_figs
          

4.3.2. The results visualization widget view

When working in Jupyter Notebook, referencing a selection of series will cause a data visualization widget to be generated that allows you to look at results (basesdf vs latestdf) for the selected variables as tables or charts, as levels, as growth rates and as percent differences from baseline.

mpak['PAKNYGDPMKTPCN PAKNYGDPMKTPKN PAKGGEXPTOTLCN PAKGGREVTOTLCN PAKNECONGOVTKN']